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The surface of conservative coupled sandpiles in the self-organized cooperative critical state is 
found to exhibit intermittency in both time and space. The spatiotemporal intermittent structure 
is also found to be a multifractal. The probability density of sand grain heights on the surface is an 
asymptotic power law but with an exponential cut-off. The power spectra of the time series of sand 
grain heights show a diversity of low-frequency components over different sites on the surface and 
also over different ensemble samples. This means that the long-term prediction according to the 
nearby observations and the history experiences is very difficult in the world of coupled sandpiles. 
Unlike the low-frequency spectra, the high-frequency spectra seem to obey a universal f~^ law. 
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In 1980s, Bak, Tang and Wiesenfeld proposed the self- 
organized criticality (SOC) theory to explain the ubiq¬ 
uitous power laws in nature [lH3[. This theory is sim¬ 
ply illustrated by a single sandpile model (BTW model). 
Recently, we find a new kind of self-organized critical 
state by coupling two BTW sandpiles Q. We call it self- 
organized cooperative criticality (SOCC) because each 
sandpile in this model is not critical but they cooperate 
to be critical. In the SOCC state, we find that the distri¬ 
butions of the avalanche size and the lifetime exhibit the 
same finite-size scalings as those in the BTW sandpile. 
This conclusion has been strictly proved to be true Q. 
However, there is a remarkable difference between the 
two kinds of criticality. In the BTW model, after each 
avalanche the negative value of local height gradient (for 
simplicity, this term is replaced by the “local height” in 
the following text) at any site will not exceed the critical 
threshold Zc because once the local height exceeds Zc an 
avalanche will be trigged again. However, in one of the 
coupled sandpiles the local height at each site may exceed 
Zc and fluctuates irregularly in space and time because 
the trigger of avalanches is defined by the two sandpiles. 
The difference between the BTW model and our model 
is significance if we turn our eyes upon the natural world. 
For example, although the BTW model can explain the 
scaling behavior of earthquakes (avalanches of sandpile) 
successfully B 0, it can not explain the very irregu¬ 
lar configurations of the earth’s surface (configuration of 
sandpile after avalanches). The irregular fluctuations in 
space and time is called the spatiotemporal intermittency 
here. In this paper, we will exhibit the quantitative fea¬ 
tures of spatiotemporal intermittency in the conservative 
coupled sandpiles. The results are expected to be used 
as fingerprints to search the SOCC state in real systems. 

The conservative coupled sandpiles are composed of 
two BTW sandpiles which are coupled in such a way that, 
if and only i/both local heights at the same coordinate of 
the coupled sandpiles (denoted by zi and Z 2 ) are greater 
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than the critical threshold Zc (equaling 3 in the following 
simulations), the sand grains then tumble down as 

z^{x,y) ^ z^{x,y) - 4, fori = 1,2 (1) 

and the tumbling down grains are transported to the 
nearest-neighbor grids, 

Zi{x±l,y) ^ Zi{x±l,y) + 1, (2) 

Zi{x,y±l) ^ Zi{x,y±l) + 1. (3) 

Note that the coordinate of the upper left corner of the 
sandpile is set to (1,1) and the coordinate of the lower 
right corner is set to {N, N) in this paper. The coor¬ 
dinates range as 1 < x < iV and 1 < y < N. The 
coupled avalanche dynamics are inspired by the coupling 
phenomena in many real systems. For example, in our 
daily life you may have considered to migrate to a big¬ 
ger city for a better job. Migration to a bigger city is 
a movement along the population gradient. Here, some¬ 
one’s final choice is determined by the equilibrium be¬ 
tween the population gradient and other gradients, such 
as the job gradients, the education gradients and even 
shopping gradients. In other words, the dynamics of mi¬ 
gration is coupled with other dynamics. In fact, similar 
coupling also exists in the natural systems. A well-known 
example is the cross effect in the transport arising in a 
mixture if both the concentrations and the temperature 
are non-uniform over the system Q. 

Starting from two sandpiles with the same size of N 
and zero local height in the whole sandpiles, the critical 
states are built by adding one sand grain randomly to 
each sandpile simultaneously at every turn. After adding, 
the two coupled sandpiles are allowed to tumble down (if 
needed) according to Eqs. ([T]), (2)) and This model 
is conservative except in the boundary where tumbled 
sand grains can flow out of the sandpile. Once the sand¬ 
piles become calm, another sand grain is added to each 
sandpile. Because of the symmetry of the two coupled 
sandpiles, we just show the statistics of one of them in 
the following analysis. 

With the increase of number of adding grains n, a crit¬ 
ical state will be built in the coupled sandpiles where 
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FIG. 1. Time evolutions of the coarse-grained outflow (top) 
and the local height averaged over the whole sandpile (bot¬ 
tom) for different sandpile sizes (circles for = 20, squares 
for = 40 and triangles for N = 60). 

the outflow will fluctuate around a constant value and 
the distributions of avalanche size and lifetime will fol¬ 
low the same power laws as those in the classic BTW 
sandpile. However, in the conservative coupled sandpile 
the local height averaged over the whole sandpile z will 
continue to increase with n in the critical state while in 
the BTW sandpile it fluctuates around a constant Si- 
It means that the surface of coupled sandpiles shaped by 
the avalanches is growing with the elapse of time. One 
can image that in larger sandpiles more adding grains are 
needed than in smaller sandpiles to form a surface with 
a defined value of z. Thus, the number of adding grains 
is not convenient to compare the degrees of growth with 
different sandpile sizes. We define the growth time tp to 
solve this problem 


where N is the sandpile size and n is the number of 
adding grains. By using tp, one can see that the curves 
of outflow measured with different sandpile sizes almost 
overlap each other (see the top plane of Fig. [J). The 
outflow begins to saturate at tp ss 4 where the sand¬ 
piles begin to be critical. By using tp, the curves of the 
averaged local height z also overlap each other (see the 
bottom plane of fig.[I]). It means that at the same growth 
time tp the degree of growth will be the same whatever 
the sandpile sizes are. In the following sections, we will 
show that tp can also be used to scale other variables into 
universal functions. 

Intermittency in space. If the local height at a co¬ 
ordinate in one of coupled sandpiles is smaller than Zc, 
the local height at the same coordinate in another sand¬ 
pile could be accumulated to be greater than Zc- This will 
lead to the spatial intermittency on the surface of cou- 


FIG. 2. A snapshot of the local heights z on a SOGC sandpile 
with A^ = 60 (blue bars for z > 20, green bars for 3 < z < 20 
and red bars for z < 3 ). The growth time tp = 50. 

pled sandpiles (see Fig. [5] for example). The quantitative 
features of spatial intermittency in the SOCC states, in¬ 
cluding the probability density functions of local heights 
and the sequence of mass exponent, are analyzed as fol¬ 
lows. 

The probability density functions of local heights are 
calculated for different sandpile sizes and different growth 
times. In order to improve statistics, we generate many 
pairs of coupled sandpiles for each sandpile size (500 
pairs for N = 20, 126 pairs for = 40 and 56 pairs 
for N = 60). We find that the probability density func¬ 
tions P{z; tp) behave as a universal function whatever the 
sandpile sizes are (see Fig. El). The universal function is 

P{z;tp) = z-^f (5) 

where a w 0.4 and (3 l.I. When is very large, 
P{z;tp) will approach to a power law but with an expo¬ 
nential cut-off, 

P(z; tp) ~ z-^ exp ^ (6) 

where ci Ri 1.5. 

The success of the BTW sandpile model may give us a 
misleading impression that complex systems always fol¬ 
low asymptotic power-law distributions such as the distri¬ 
butions of earthquake sizes (the Gutenberg-Richter law 
0), the incomes (the Pareto’s law M) and the word 
frequency (the Zipf’s law [Hj). However, there are many 
other complex systems which do not follow the asymp¬ 
totic power laws [l^. In fact, the asymptotic power law 
with an exponential cut-off as one kind of the non-power- 
law distributions are observed in many systems. Exam¬ 
ples include the turbulent wind fluctuations in the atmo¬ 
spheric boundary layer [mil, the complex biological 
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FIG. 3. Probability of density functions of local heights 
P{z\ tp) for different sandpile sizes (black points for N — 20, 
red points for Ai = 40 and blue points for N = 60) and growth 
times (circles for tp = 10, squares for tp = 15 and triangles 
for tp = 20) in the SOCC state. The slope of line is about 
-1.5. 


[l^ and social [13 networks. Equation ([S]) shows that 
the SOCC could be an alternative mechanism of this kind 
of non-power-law distributions. 

The sequence of mass exponents is a useful tool to de¬ 
tect and describe the multifractal structure. Let us cover 
the whole sandpile with boxes of length S and calculate 
the probability of local height in the i-th box 


fJ-i = 


E m 

i=l^i 


( 7 ) 


where m is the least number of boxes needed to cover the 
sandpile. If the sum of g-th moments of jii has the form 


Z{q,5) = ynl 


2=1 


J-r(g) 


( 8 ) 


when (5—^0 and if the power index r(g) is a non-linear 
function of q, the surface of sandpile is then a multifractal 
[T^ . Equation ([5]) is tested to be true for different growth 
times in the SOCC state (see the left panel of Eig. H. The 
mass sequence r(g) is also found to be a non-linear func¬ 
tion of q and can be fitted well with the mass exponents 
of 2D binomial multiplicative cascade process |1^ . 


-(?) = 


In (c 




(P) 


In 2 


(9) 


where a = 0.3, b = 0.3, c = 0.23 and a + b + c + d= 1 
(see the right panel of Eig. n. The results show that 
the surface of coupled sandpiles in the SOCC state is a 
multifractal and its self-similarity feature will not change 
with tp. 

Intermittency in time. The time series of local 
height at different sites are shown in Eig. Many in¬ 
teresting features of coupled sandpiles are shown in this 
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FIG. 4. (Left) The sum of g-th order moments of probability 
of local height Z{q,S) as a function of box length S. Lines 
show that the data with different values of q can be well htted 
by the power law when S is not very large. From the bottom 
line to the top line, the order q = —15, —10, —5, 5,10, 15. In 
this plot, the side length of the square sandpile is set to 1 
and the sandpile size N = 60. (Right) The sequence of mass 
exponents r(g) as a function of moment order q. The line is 
calculated from Eq. (|9]). In both panels, circles, squares and 
triangles denote the results for tp = 10, tp = 15 and tp = 20 
respectively. 


figure. First, the coupled sandpiles have the complicate 
intermittent patterns in time even the sandpile size N is 
not very large. The time evolutions of local height with 
a size of only TV = 10 show an variety of patterns (see 
Fig. [5]). Second, adjacent sites could have dramatically 
different patterns of time evolution. For example, one 
can see that the dramatic difference between the time 
evolutions at adjacent sites (6, 5) and (6,4) (see the mid¬ 
dle and bottom planes of Fig. 0. This means that we 
can not make a good prediction of nearby sites according 
to the local experiences in the coupled sandpiles. Third, 
sites that are separated further from each other may have 
a similar pattern of time evolution. For example, in the 
top and middle planes of Fig. [Sj one can see that both 
time series at the sites (9,10) and (6, 5) have a long-term 
increasing trend first and then saturate at tp ^ 600. We 
are not sure whether this similarity is related to some 
physical connection or just a coincidence. 

The power spectra S{f) of the time series of local 
height are analyzed at every site in the SOCC states. In 
the analysis, the sampling interval is set to Atp = 1/A^^. 
Results show that the high-frequency spectra are similar 
for all sites and behave as a power law, 


Sif) ~ / 


-« 


( 10 ) 


where ^ « —2 (see the left panel of Fig. 0. This seems 
to be true for different sandpile sizes. However, the low- 
frequency spectra are very different. This means that 
the long-term variations in the time series of height are 
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FIG. 5. Time series of local height at sites (9,10) (top), (6, 5) 
(middle) and (6,4) (bottom). The sandpile size N = 10. 
Note that the time series at the same site will change with 
the ensemble sample of coupled sandpiles. The time series 
shown here are from a possible ensemble sample. 



FIG. 6. (Left) Bin-averaged spectra at all sites on the SOCG 
sandpile. Green lines are for the sandpile size = 20 and 
blue lines are for N = 10. (Right) Bin-averaged spectra at 
the site (20, 20) on the SOCC sandpile with a sandpile size of 
N — 40. Spectra of 50 ensemble samples are shown in this 
plot. 


diverse over the sandpile which is a reflection of the spa¬ 
tial intermittency. We generate many pairs of coupled 
sandpiles with a defined sandpile size and then analyze 
the spectra at a defined coordinate. It is found that the 
high-frequency spectra can also be described by a power 
law of Eq. (fTOl) for each ensemble sample. However, the 
low-frequency is also diverse over different ensemble sam¬ 
ples (see the right panel of Fig. IH]). This means that the 
randomness will severely affect the long-term variations 
in the coupled sandpiles. Thus, it will be very difficult 
to make a long-term prediction according to the history 
experience in the coupled sandpiles. 

In conclusion, we have analyzed many quantitative fea¬ 
tures of spatiotemporal intermittency on the growing sur¬ 
face of coupled sandpiles in the SOCC state including the 
probability density functions of local height, the sequence 
of mass exponent and the spectra. It is found that these 
features exhibit some universal regularities, although the 
surface of coupled sandpiles seems to be very complicated 
in space and time. These findings may help us to diag¬ 
nose SOCC states in real systems. 
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